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Abstract 

A systematic diagrammatic expansion for Gutzwiller-wave function (DE-GWF) is used for the 
description of superconducting state in the Hubbard model. The method ahows for an analysis of 
the superconducting-gap evolution between the weak and the strong-coupling limits. The super- 
conducting gap at the Fermi surface resembles d-wave only near and at the optimal doping and 
the corrections to this state are determined; they are coming from the longer range nature of the 
resulting pairing. The superconducting state appears only for U/t > 3 and for the doping 6 < 0.32; 
the latter agrees very well with experimental results for the cuprates. The spontaneous breakdown 
of the Fermi-surface coherence (formation of the Fermi pockets) in the antinodal direction is also 
discussed. General effective Hamiltonian with pairing is defined in a fully self-consistent manner. 

PACS numbers: 71.27.+a, 74.20.-z, 74.20.Rp 



1 



High-temperature superconductivity is very often discussed starting either from the Hub- 
bard model l| or from its projected version in the strong-correlation limit - the t-J model 
These models incorporate, in the simplest manner, the strongly correlated nature of 
the 3(i electrons due to Cu ions in CUO2 planes. While within the t-J model one intro- 

n 

duces real-space operators for an antiferromagnetic coupling explicitly [3[, the situation is 

less obvious in the former case unless one introduces antiferromagnetic spin-fluctuations as 

fl 

mediating the pairing jj], which can be analyzed reliably for low or moderate values of the 
Hubbard interaction U . In general, we require a method interpolating between the original 
version and its projected t-J version to test the universality of this nonstandard, real-space 
pairing concept, as well as to interpolate between the weak and the strong correlation limits. 
Only then these nonstandard ideas will be placed on a firm conceptual ground. Also, the 
pairing evolution with the increasing interaction strength is particularly interesting in view 
of the circumstance that the recently found iron-pnictide superconductors are regarded as 
moderately correlated systems 5[. The Variational Monte Carlo (VMC) method is one of 
such methods limited however to small-size systems, containing typically 10 x 10 lattice 
sites in the two-dimensional situation js, 7|. Comparable accuracy (and limitations) have the 
works generalizing the density matrix renormalization group approach applied to t-J model 
jsj. An extensive analysis has also been carried out within the 2x2 cluster dynamical 
mean-field theory joj. 

We propose the extension of a recently devised systematic expansion based on the Gutz- 
willer Wave Function (GWF) method for the paramagnetic state 10|], which we apply to 
determine stable superconducting state and its systematic evolution with the increasing U 
and doping 5. In particular, we determine the doping dependence of the superconducting 
gap and also show that the anisotropic (k-dependent) gap at the Fermi surface has rf-wave 
character only around and at the optimal doping. The universal character of the disappear- 
ance of the superconducting gap at the doping 5 ^ 1/3 is established, in agreement with the 
recent sophisticated renormalized mean- field theory results [ll| and experiment. We also 
calculate the longer range components of the gap parameter Ai j in real space (i.e., beyond 
the nearest-neighbor contribution A^ij^) and note their systematic evolution, with increase 
of the doping. Finally, we note a systematic difference with the doping dependence of the 
Fermi surface topology in the antinodal direction with respect to that in the nodal direction 
and term it "breakdown of the Fermi surface coherence" . 
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The GWF method has been tested in one spatial dimension 12|, where it removes the spu 



rious metal-insulator transition present in the Gutzwiller approximation (GA) and compares 



favorably with the exact Lieb-Wu solution 13|. In this letter, we apply the Diagrammatic 



Expansion of GWF (termed DE-GWF) approach to the case of a two-dimensional square 
lattice with two hopping integrals (— t) and t' = 0.25 1 and concentrate on the supercon- 
ducting (SC) solution. The main features of the DE-GWF method for the paramagnetic 



(normal) state have been provided in Ref. 
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Here we give a brief overview of the main 



steps and subsequently generalize the approach to the description of the SC case. We start 
from the Hubbard Hamiltonian in the form 

H = Ho + U^di, Ho = ^ tijcl^c.,^ , di = rii.t'^i,;, (1) 

where i = (11,12) is the two-dimensional site- index, tjj = —t is the hopping integral for 
nearest neighbors (n.n.) and t' for the second neighbors, whereas a =1, f is the spin quantum 
number. The Gutzwiller wave function [l^ for the correlated state has the form 

|M/G)=P|^0) = ni^i|'^0)' (2) 

where |\E'o) is a single-particle Slater-product state to be defined later. In this work we use 
the local projector as defined by 



A = $^Ar|r)H(r|, (3) 

r 

P^ = l + xd?^. (4) 



In Eq. ([3]) a general form of the projector has been presented, with the variational parameters 
Ar € {A0, All, Aij^, A(i} describing occupancies of the four possible local states {|r)i} = 
1 t)i) I I ti)i}- III Eq- ©) a particular form of the projector is assumed, with 
the Hartree-Fock operators defined by df^^ = "^^i^^^^i^^^ and = Ui^cr — no, with no = 
(^o|^i,cr|^o)- Such form of simplifies essentially the calculations by eliminating the 



'Hartree bubbles'' 
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15|. In this manner, three out of the four parameters {Ar} are fixed. 
The remaining of them can be either used as a variational parameter or expressed via x 
and then x is selected as a variational parameter, which turns out to be a computationally 
efficient choice. Next, having selected the form of the projector and the variational parameter 
X, we can calculate all the needed expectation values: the norm, (\E'g|^g)5 the double 



occupancy (\l/GMi|^G), and the hopping term (^g|c; o-^j o-I^g)- Details of the calculations 



can be found in Ref. 



lOt the summary has been provided for the sake of completeness 



in the Supplemental Material |l6j. Here we concentrate on discussing the new features 
appearing in the case of superconducting phase. First, apart from the "normal" |\Efo) lines, 
as represented by Pi y = PC^, = (cj^Cj, ^)o — SiyriQ, we now also have the "superconducting" 
lines S'lj/ = (cj^Cj, ^)o. Note that we will only consider superconducting orders without 
a local pairing, i.e., with ^i^i = 0. Second, since the correlated number of particles, uq = 
(^i,<T)G and its non-correlated correspondent no may differ in the superconducting phase, the 
minimization procedure is different. Namely we minimize the generalized grand-canonical 
potential J-" = {H)q — 2^g^gL instead of minimizing the ground state energy Eg = {H)q 
Third, the minimization procedure leads to an effective single-particle Hamiltonian which, 
in the present situation, contains also the superconducting pairing part, as defined by 



where f =i and l =t. The above form of the Hamiltonian allows for definition of renormal- 
ized masses which will be discussed below. Finally, an interesting additional quantity in the 
present case is also the correlated gap, defined by Ac = {cl-^cj^)G (see Supplemental Material 
for the explicit analytical expression) for nearest neighbors i, j. The solution procedure for 
the superconducting case, together with the resulting self-consistency loop, is also presented 
there. 

We underline here that the analysis of the superconducting case is much more involved 
computationally than that for the paramagnetic case, because the number of diagrams grows 
substantially (e.g. in the 5-th order there can be 1000 times more "superconducting" dia- 
grams than "normal" diagrams). To be able to perform calculations for the superconducting 
phase in the 5-th order of the expansion, an additional improvement had to be made. The 
improvement relies on counting the equivalent contributions when performing a summation 
and thus obtaining the analytic forms of the diagrammatic sums (These forms can have 
~ 6 X 10^ components). Such scheme requires the use of the hash table data structure to 
store information concerning different contributions. 

Now, we discuss the physical results coming from the implementation of our DE-GWF 



method. Except when stated otherwise, the presented results were obtained in the 5-th 
order of the expansion and for the values of parameters: t' = t/4, and for t = 1 (which 
sets the energy scale). We also take into account the |\l/o) lines Pij = -Po,(i-j) = Pxy (with 



X = ii - ji, Y = 12-32) fulfilling X'^ + Y'^ < 10 (and the same condition for ^ij, tfj, A??^ 
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FIG. 1. (color online), a) Phase diagram containing paramagnetic (PM) and superconducting 
(SC) phases on plane interaction strength (U) - doping (5) for selected values of t'. The reentrant 
behavior is associated with the dome-hke SC appearance on phase diagram as a function of 6 (cf. 
the U = 6 curve in (c)). b) Condensation energy (in Kelvins, assuming t = 0.35 eV) as a function 
of doping for selected U values. The inset shows the gain in kinetic (Aii^kin) and potential (AE'pot) 
energy with respect to the paramagnetic state for U = 10. c) Correlated gap vs 5, with well defined 
upper critical concentration 6 ~ 1/3. d) Correlated gap vs U for selected dopings. 



In the panel composing Fig. [T] we display the obtained principal ground-state characteris- 
tics of the superconducting phase, namely, the phase diagram (a), the condensation energy 



(b), defined as AEq = e\^'^^^ —Eq^\ the correlated gap vs doping 6 = 1 —2nG (c), and 



n(SC) 
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Ag vs. U (d). The superconducting region in a) expands with the increasing t' towards the 
higher-doping values, whereas the onset of the paired-state appearance for f/ > 8 is practi- 
cally independent of U. The practical advantage of our method is that the phase diagram 
in Fig. [T^ can be computed in time scale of a few days on a modern PC. The condensation 
energy shown in Fig. [TJj is in Kelvins (assuming t = 0.35 eV), to show that this quantity is 
in accordance with the experimental values for the cuprate high-Tc superconductors. The 
inset in Fig. [T|d shows the corresponding changes in the kinetic (AE'kin) and the potential 
(AE'pot) energies counted with respect to the paramagnetic state; the superconductivity at 
[/ = 10 is due to a gain in the potential energy. This is not always the for higher U 

values (e.g. ai U = 14) superconductivity is driven by a gain in the potential energy (close 
to the half filling), both energies (at dopings 0.04 < 6 < 0.09), and the kinetic energy (at 
dopings 6 > 0.09). In Fig. [Th the correlated gap shows a dome-like behavior for U > 10. One 
of the fundamental features of our results (cf. Fig. [T|d) should be emphasized. Namely, the 
maximal value of the condensation energy is achieved at ?7 ~ 10, a representative value for 
a cross-over from moderately to strongly-correlated regime. This result may suggest, where 
to look for the material with the highest value of Tc. Also, the condensation energy is very 
small in the weak-correlation limit U < 6. Thus, this type of paired state should not appear 
in the weak-correlation limit. The small condensation energy for f/ < 6 means that the 
exact shape of the phase diagram in this regime depends on what we define as the onset of 
the superconducting phase. Namely, we have regarded the phase as superconducting when 
A^Q > 10~^. The superconducting state disappears for the doping 6 ^ 0.32, in very good 



agreement with both sophisticated version of the renormalized mean-field theory Uj and 
with experimental data for virtually all single-plane cuprates. 

To visualize one of the unique features of the present approach we plot in Fig. |2] the 
effective components of the g cip clS db function of the doping. The plots of the parameters 
A;? = A^fy (with X = ii — ji and Y = i2 — 32) in Fig- show that the dominant 
component is A^q (n.n. contribution), hence the gap has mainly the (i^2_j^2-wave character. 
Ability to determine the more distant components of the gap is another principal feature 
of our approach. We also plot in Fig. Wp the effective gap in the k-space across the Fermi 
surface. It turns out that this gap deviates from the (ia:2_^2-wave form, particularly in the 



antinodal directions. Such deviation has been observed in high-Tc superconductors 17H20|. 



but we cannot claim decisively that our results reflect the physics of this phenomenon, as 
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FIG. 2. (color online), a) The effective gap parameters as a function of doping for U = lOt. 
Inset: the gap relative to the dominant A^g component, b) The effective gap in the k-space at the 
Fermi energy for a choice of doping values and U = Wt. The thin black line is a guide to the eye 
indicating momentum dependence of the standard (1^2 _y2 form. The gaps are normalized so that 
the d^2_y2 component has the value of 1 in the antinodal direction. 



experimentally the deviation is assigned to be caused by two energy scales corresponding to 
two gaps present , whereas in our case it comes from the nonzero contribution to the gap 



A?? from more distant neighbors. Note also, that the gaps A?? presented in Fig. [2] enter in 
the effective dispersion relation of the Hamiltonian ([5]), hence they can be related to, e.g., the 
photoemission experiments, whereas the correlated A^ is the gap in the thermodynamical 
sense. In the overdoped regime {6 > 0.2) the gap components are of comparable magnitude, 
what might indicate a gradual transition from real to reciprocal (k-) space pairing as the 
Fermi-liquid regime is reached. 
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FIG. 3. (color online). Mass enhancement for the paramagnetic state at the Fermi energy (as a 
function of azimuthal angle in the first quarter of the Brillouin zone) for selected doping values. 
Inset: ratio of the maximal mass enhancement Q^^^ (typically in the antinodal direction), to that 
in the nodal direction (5^'^^°"^^^ vs 6. 

It is also important to provide principal characteristics of the normal phase. In Fig. 12] we 
analyze the relative mass enhancement for the paramagnetic state, defined as Qk = /m^^. 
For both masses we use oc l/(|9kek|) calculated at the Fermi surface. For the mass 
we use the effective dispersion of the Hamiltonian ([5]), whereas for we take the bare- 
band dispersion of the Hamiltonian ([T]). It can be seen that the mass in the antinodal 
direction is much larger (even over 2 times) as compared to that in the nodal direction. 
This observation means that within the present approach we observe a tendency towards 
the pseudogap formation. Namely, the mass in the antinodal direction may be exceeding the 
critical value for the Mott localization (especially at 5 ~ 0.01), thus prohibiting a coherent 
hopping in this direction (such situation already takes place in the c direction). On the 
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other hand, the mass in the nodal direction may be substantially smaller, enabling thus 
itineracy in this direction. Therefore, the present results are in agreement with the angular 
mass dependence of high-Tc superconductors observed e.g. in ARPES experiments as Fermi 



pockets 



21 



22| . The present pseudogap-formation tendency is rooted solely in the strong 
electronic correlations; this argument contributes to the vital discussion of the possible 
pseudogap-formation mechanisms 



For completeness, in the inset we show the ratio of 
the maximal mass enhancement to the mass enhancement in the diagonal direction. This 
ratio increases with approaching half-filling and with an increase of the interaction strength 
U. 
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FIG. 4. (color online), a), b) Comparison with the Variational Monte Carlo results: doping 
dependence of the n.n. gap parameter (a) and of the condensation energy (b). The data points 
(three on each graph) are taken from Ref. The lines show our results: the solid lines represent 
a VMC-like results obtained within the present approach (for details see main text), whereas the 
dashed lines show our full results, c), d) Double occupancy vs interaction at half- filling for the 
paramagnetic (c) and superconducting (d) phases. The critical value at which do = is about 
Uc ~ 14. The numbers 0-5 label the order to which the expansion is carried out. 

Finally, we compare in Fig. H] the results of our method to the VMC results of Ref. Q, 
obtained for U = 8 and t' = 0. The "VMC-like" curves were obtained by the present DE- 
GWF method (in the 5-th order) by setting the effective parameters t?? and Af? to zero 
for all i,j's except for n.n. The remaining t^f was set as equal to —t, whereas A^q = A 



was chosen as a minimization parameter 16|. The "DE-GWF" data are the result of our 

u 

full 5-th order expansion. The VMC results of Ref. [71, and the VMC-like results of our 
method are close to each other at half filling, but qualitatively different away from it. The 
reason for this is probably that the 8x8 lattice used in the VMC calculations is too small 
to emulate the infinite-lattice size we are considering. We have also checked that increasing 
the number of |\l/o) lines {Pxy, Sxy) to those with + < 16 does not yield qualitative 
changes (in the 3-rd and 4-th orders). Also, the 3-rd, 4-th, and 5-th order results are close 
to each other, and therefore the discrepancies in the VMC and DE-GWF. results cannot be 
due to a slow convergence of our method. The right panel shows the doping dependence of 
double-occupancy probability for orders 0-5 to which the expansion is carried out. 

Summarizing, we have successfully formulated an efficient Diagrammatic Expansion 
method (DE-GWF) based on the Gutzwiller Wave Function approach, that allows for a 
systematic analysis of the superconducting-state evolution as a function of either the doping 
or the magnitude of the local interaction U. Note that the gap at the optimal doping 
reaches its highest value near U ~ 10. Also, a longer-range nature of the superconduct- 
ing correlations, as expressed via the effective gap function A??, leads to the corrections 
to the pure rf-wave type of pairing. The 6 dependence of more distant gap components 
may suggest a gradual evolution from the real-space-pairing limit [6 < 0.15), where n.n. 
component dominates, to the BCS-like limit (for 6 > 0.2), since in the latter regime the 
gap components are comparable, albeit quite small. We have also reproduced the angular 
mass dependence observed as Fermi pockets in ARPES experiments for the normal phase 



of high-Tc systems [2l|, |22| . In general, the Hubbard model with our GWF expansion leads 
to a viable superconducting ground state in moderate- and strong-correlation limits. In 
this manner, we have extended the discussion of the paired state in the Hubbard model on 
infinite lattice beyond the currently used versions of renormalized mean-field theory 
and VMC analysis A competition or coexistence with either antiferromagnetic or 

Pomeranchuk phases, as well as the extension to multi-band systems can be also formulated 
and should be investigated in the near future. Finally, a direct comparison of the results 
obtained here with those of the t-J model within the same approach is highly desirable, but 
requires a separate formulation of the corresponding diagrammatic expansion for that case. 
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